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Abstract 

FEpX is a modeling framework for computing the elastoplastic deformations of polycrystalline 
solids. Using the framework, one can simulate the mechanical behavior of aggregates of crystals, 
referred to as virtual polycrystals, over large strain deformation paths. This article presents the 
theory, the finite element formulation, and important features of the numerical implementation that 
collectively define the modeling framework. The article also provides several examples of simulating 
the elastoplastic behavior of poly crystalline solids to illustrate possible applications of the framework. 
There is an associated finite element code, also referred to as FEpX, that is based on the framework 
presented here and was used to perform the simulations presented in the examples. The article serves 
as a citable reference for the modeling framework for users of that code. Specific information about 
the formats of the input and output data, the code architecture, and the code archive are contained 
in other documents. 


*http://anisotropy.mae.cornell.edu/dplab/ 
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1 Introduction 


1.1 Purpose 

The purpose of this article is to lay out a complete system of equations for modeling the anisotropic, 
elasto-viscoplastic response of polycrystalline solids comprised of aggregates of grains and to summarize 
a finite element formulation that may be employed to compute the motion and stress in polycrystals 
governed by the system of equations under imposed loadings. The governing equations together with 
associated solution methodologies define a modeling framework, referred to as FEpX, that is focused 
at a physical length scale of an ensemble of grains. There is an associated finite element code, also 
named FEpX, that follows the framework. A major motivation for archiving this article is to provide 
a thorough and accessible reference that researchers who utilize the code can readily cite. However, 
the article stands independently in providing a complete summary of a crystal-scale model for the 
elasto-viscoplastic response of polycrystalline aggregates and a finite element formulation that enables 
solving the model equations over motions that entail large deformations. 

The content provided here regarding the governing equations and finite element framework draws 
primarily from the following published articles: PEIEIII]. The present article is not intended to serve 
as a primer for computational crystal plasticity, so background knowledge of solid mechanics, including 
crystal plasticity, and nonlinear finite element methods is assumed. Rather, it strives to encapsulate 
the full set of equations, assumptions, and solution approximations necessary to document simulation 
results with sufficient detail to facilitate those results being reproduced by others. 

1.2 Scope 

The scope of this article is limited to the theory and methods that define the FEpX framework, plus 
a general overview of the data flow within the framework and the interfaces with tools to instantiate 
virtual polycrystals and to visualize simulation results. Consequently, there are sections of the article 
devoted to these topics, as listed in the Table of Contents. Also provided are representative examples to 
illustrate application of the derivative finite element code to modeling of single and dual phase metallic 
alloys. No detailed information is included on the specific formatting used for problem definition, code 
execution, or exported simulation results. That information is contained in separate documentation 
associated with the code itself. 

1.3 Complementary modeling tools 

The role of FEpX in the modeling of polycrystals is to solve the boundary value problem associ¬ 
ated with the elastoplastic response of a polycrystalline solid arising from applied mechanical loading. 
Separate tools are needed to instantiate a virtual polycrystal and to discretize it with finite elements. 
FEpX accepts the finite element mesh generated by custom MATLAB scripts (available in the 
OdfPf package) and by the Neper program [5]. Separate packages for visualizations also are needed. 
Export scripts are available for writing files that can be imported by Paraviewjn], and Visits . 
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2 Capabilities of FEpX 


The FEpX framework is a combination of a theoretical construct developed by a world-wide commu¬ 
nity of researchers and a numerical solution methodology (finite element formulation) developed by the 
DPLab members over the past than two decades. The governing equations and solution methodology 
become tightly intwined through the choice of interpolation functions for the motion and the weighted 
residuals for equilibrium. Thus, we refer to the combination as a framework and do not make a strong 
distinction between the framework and the derivative code, calling both FEpX. 

2.1 What FEpX can do. 

FEpX is useful for simulating the mechanical behavior of polycrystalline solids at the level of aggregates 
of grains. The aggregates may be comprised of grains of a single phase or of multiple phases. Grains 
are discretized with finite elements so any sub-volume of an element is a sub-volume of an individual 
crystal. The local behaviors associated with the material with any element correspondingly are those 
of a crystal. In particular, the behaviors include: 

• nonlinear kinematics capable of handling motions with both large strains and large rotations; 

• anisotropic elasticity based on cubic or hexagonal crystal symmetry; 

• anisotropic plasticity based on rate-dependent slip on a restricted number of systems for cubic 
or hexagonal symmetry; 

• evolution of state variables for crystal lattice orientation and slip system strengths; 

To accommodate these behaviors the finite element formulation has incorporated a number a numerical 
features, such as; 

• higher-order, isoparametric elements with quadrature for integrating over the volume; 

• implicit update of the stress in integrations over time; 

• monotonic and cyclic loading under quasi-static conditions; 

Depending on the goals of the simulation, aggregates might number in grains from only a few to 
tens of thousands (or more). The grains can be discretized at a level appropriate for the intent of 
a simulation. The number of grains together with the level of discretization within grains set the 
computational burden for a simulation. To accommodate combinations with high burden, the FEpX 
code employs scalable parallel methods and executes on clusters. FEpX has been developed to use 
meshes constructed by instantiation tools for virtual polycrystals and to output data in an archivable 
format for subsequent use with visualization tools or other interpretation tools. 

With these capabilities FEpX is well-suited, for example, to model the mechanical behavior of 
polycrystals that exhibit inhomogeneous deformations within and among the crystals, to investigate 
the heterogeneity of stress within a polycrystal, or to examine the roles of neighbors on the behaviors of 
individual grains. When teamed with appropriate instantiation methods, EEpX can be used effectively 
to model yielding and flow of alloys with complicated phase/grain topologies and morphologies. 
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2.2 What FEpX cannot do. 

The FEpX framework does not encompass many aspects of the behaviors observed in real materials. 
Some of its principal limitations include: 

• plastic flow occurs by slip - no other mechanisms, such as twinning and creep, are modeled; 

• deformations are ductile - no fracture models are included; 

• loading is quasi-static - no inertial effects are modeled; 

• loading is mechanical (isothermal) - coupling with heat transfer (or other physical processes) is 
not considered; 

• boundary conditions are simple - neither friction models nor changing contact conditions are 
included. 

With these limitations, FEpX is not well-suited for modeling applications with complex loading con¬ 
ditions, such as many metal forming and joining processes, or for modeling applications involving 
fragmentation failure of a dynamically loaded body. 

2.3 What FEpX has been used to study. 

A variety of interesting problems arise at physical length scales in which a sample volume encompasses 
an aggregate of grains. We typically think of an aggregate containing on the order of 10^ — 10® grains, 
but the simulation framework embodied in FEpX is appropriate for single-grain or multi-grain samples 
(1 — 10^ grains), as well. Some of the applications of FEpX published in the open literature are listed 
below. Users of FEpX are encouraged to examine articles in areas of interest to obtain information 
beyond the scope of this article resulting from the collective experiences of others in applying FEpX. 

• Grain interactions with attention focused on bulk texture evolution. Articles published in this 
area are: m El uni nn da nail!. 

• Deformation heterogeneity within the grains comprising an aggregate with focus on intra-grain 
misorientation distributions. Articles published in this area are: da da da ds]. 

• Inter- and intra-grain stress/elastic strain distributions, especially including comparisons to neu¬ 
tron and x-ray diffraction experiments. Articles published in this area are: daEanaEaEaiMj. 

• The elasto-plastic transition occurring during loading of polycrystalline solids, with focus on the 
redirection of stress at the grain level. Articles published in this area are: [HEaEaEa- 

• Cyclic loading with interest in evolution of stress and its implications for fatigue failure. Articles 
published in this area are: p8l[29] . 

• Evolution of dislocation density and associated peak broadening. Articles published in this area 
are: [HOlEI]. 

• Virtual polycrystal instantiation issues, including sensitivity of the stress and deformation to 
discretization. Articles published in this area are: [33 ESI E]. 
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Much of the original work in utilizing polycrystal plasticity constitutive models within finite element 
simulations was focused on bulk texture evolution in macroscopic scale deformation processes, such as 
metal forming operations (rolling, extrusion, and sheet forming) and geological processes (particnlarly 
mantle convection). In such cases, the relative sizes of finite elements and grains were reversed in 
comparison to those of the intended applications of the FEpX framework described in this article. 
Conseqnently, the mechanical properties within an element were derived from an average over an 
ensemble of crystals once an averaging hypothesis {e.g. isostress or isostrain) was imposed. Examples 
of this type of application are: [Ml |35l ESI EZl ESI ESI HQl |2T]- While relevant from a historical 
perspective in the development of FEpX, the FEpX framework described here does not include 
evaluating properties within an element on the basis an average over a population of grains. Rather, 
properties within an element are those of a single orientation, consistent with an element spanning only 
a part of any given grain. 


3 Nomenclature 


3.1 Variables used in theoretical description 

13 - determinant of the elastic stretch, 

7 “ - shearing rate of the a slip system 

u - surface normal vector 

L - body force vector 

TT - mean stress (tr(<T)/3) 

(j) - rotation angle associated with r 

a - Cauchy stress 

T - Kirchhoff stress 

r“ ~ resolved shear stress on the a slip system 
u - spin vector associated with v 

X ~ mapping function of motion 

B - domain of the polycrystal 

dB - surface of the polycrystal 

C - elasticity (stiffness) tensor 

d - deformation rate (symmetric part of /) 

dP - plastic deformation rate (symmetric part of /^) 
e® - elastic strain 

f - deformation gradient 

- elastic part of the deformation gradient 

f* - rotational part of the deformation gradient associated with the lattice rotation 
fP - plastic part of the deformation gradient 

51 “ - strength of the ol slip system 

/ - velocity gradient 

IP ~ plastic velocity gradient 

m“ - normal to the slip plane for the a slip system 

M. - plasticity (stiffness) tensor 

n - axis vector associated with r 

p" - symmetric part of the Schmid tensor for the ol slip system 

( 7 “ - skew part of the Schmid tensor for the a slip system 

q - quaternion representation of lattice orientation (go, ^ 

r - Rodrigues vector for the orientation of the crystallographic lattice 

r* - rotational part of the deformation gradient associated with the lattice rotation (= P) 

R - rotation operator corresponding to r 

s" - slip direction for the ol slip system 

t ~ traction vector 

t - imposed traction vector on the surface 

At - time step 

V ~ velocity vector of a point in the current configuration 

V - imposed velocity vector on the surface 

V® - elastic stretch 

w - spin (skew part of /) 

- plastic spin (skew part of P) 

X - position vector of a point in the current configuration 

X - position vector of a point in the reference configuration 
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3.2 Parameters appearing in the constitutive models 

7 o “ fixed-state strain rate scaling coefficient 

'js - saturation strength strain rate scaling coefficient 

K - elastic bulk modulus 

Cij - components of the elastic stiffness 

go - initial slip system strength 

gi - reference value of saturation strength 

ho - strength hardening rate coefficient 

m ~ hxed-state strain rate sensitivity 

m' - saturation strength rate scaling exponent 

n' - power on modihed Voce hardening term 


3.3 Variables used in implementation description 





B 

c 


{d} 

{en 



m 



R 


u 




{x} 


X 


- matrix trace operator 

- matrix form of the weights 

- nodal point weights vector 

~ vector matrix form of Cauchy stress 

- vector matrix form of Kirchhoff stress 

- spatial derivatives of the interpolation functions, [N] 

- matrix form of the elastic stiffness 

- vector matrix form of deformation rate 

- vector matrix form of elastic strain 

- elemental and global surface traction and body force matrices 

- elemental and global initial elastic strain matrices 

- elemental and global spin correction stiffness matrices 

- vector matrix form of the spin correction and initial elastic strain terms 

- elemental and global deviatoric stiffness matrices 

- elemental and global volumetric stiffness matrices 

- matrix form of the plastic stiffness 

- interpolation functions for interpolation of the velocity distribution 

- vector matrix form of the symmetric part of the Schmid tensor 

- equilibrium weighted residual 

- elemental residual vector 

- matrix form of the plastic spin 

- matrix form of the velocity 

- nodal point velocity vector 
~ matrix form of the position 

- coefficient matrix of factors necessary to deliver correct inner product from matrix multiplicat 
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4 Governing Equations 

The schematic diagram shown in Figure[^is a useful depiction of the intended application of the FEpX. 
The body is a polycrystal in which the full volume is subdivided into the individual crystals such that 
the entire volume is completely filled. The inter-crystal boundaries are cohesive so that there is neither 
sliding at the crystal boundaries nor separation between crystals. Every crystal is discretized with 
one or more finite elements. The material properties are evaluated at quadrature points within the 
elements on the basis of single crystal behavior, which is intimately tied to the crystal structure and 
its orientation with respect to the body at large. External load are applied to the polycrystal along 
with the appropriate kinematic constraints. 



Polycrystal 


k. 


Single 

Crystal 



Balance Laws 


Finite Slip 

Element Modes 

- 

Constitutive Model 


Eigure 1: Schematic diagram of a polycrystal subjected to mechanical loading and the discretization 
of individual crystals with finite elements. 


4.1 Kinematics and balance laws 

The motion of the virtual polycrystal is assumed to be a smooth mapping of all points within the 
domain, B. The domain in this context refers to the union of all crystal volumes that collectively 
define the virtual polycrystal. The internal crystal interfaces are smooth surfaces that remain coherent 
throughout the motion. Under the motion, the current coordinates may be written as a function of 
reference coordinates for every point in the domain; 

X = x{X) 

The deformation gradient is defined locally as: 

r _ dx 

^ 1)X 

While the mapping, y, is smooth, the deformation gradient will be so only within the interior of ele¬ 
ments. Discontinuities can arise across element boundaries, whether the boundaries lie within crystals 
or on the interface between crystals. The time-rate-change of the deformation gradient is given by: 


f = l-f 


(3) 
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where the velocity gradient, /, is computed from the velocity field, v(x), as: 

dv 


I = 


dx 


( 4 ) 


Again, while the velocity is smooth everywhere, its gradient may have discontinuities across element 
boundaries. The velocity gradient is decomposed into its symmetric and skew parts: 


I = d + w 


(5) 


where d is the deformation rate (symmetric part) and w is the spin (skew part). 

The motion of the polycrystal is driven by the stress. The Cauchy stress, cr{x), is a field variable 
defined over the polycrystal domain. Under the loading assumptions, inertia in the balance of linear 
momentum is neglected, giving static equilibrium in the local form as: 


diver"' +1 = 0 


( 6 ) 


where l is the body force vector. Body forces are neglected in the current implementation of FEpX. 
This equation applies to the interior of the crystals. Across crystal interfaces, continuity of the traction 
is needed. Applying the Cauchy formula, this condition may be written for two contacting crystals, i 
and j, as: 

t{xy = t{xy (7) 

where the tractions are related to the stress by means of the Cauchy formula: 

t = I'ix) • cr (8) 

The Cauchy stress may be split into deviatoric and spherical parts: 

a = a' — ttI (9) 

which is central to the material response as only the deviatoric part drives plastic flow. 

4.2 Constitutive equations 

The material behavior is quantified with a set of constitutive equations, here written at the level 
of the single crystal. The behavior includes both elastic (recoverable strains upon removal of the 
stress) and plastic (non-recoverable strains upon removal of the stress) responses. These can occur 
concurrently which requires coupling of the motions for the two via a kinematic decomposition. The 
decomposition is not strictly derivable from the mapping, but rather involve assumptions regarding 
the behavior. Consequently, it is part of constitutive model. The elastic response is limited to linear 
behavior following Hooke’s law for anisotropic behavior. The plastic response is nonlinear and rate- 
dependent (viscoplastic). It is assumed to be isochoric and independent of the mean stress. The set 
of equations are summarized in the following subsections, starting with the kinematic decomposition. 
The equations for the elastic and plastic responses follow discussion of the kinematic decomposition 
and are broken into two parts: fixed state relations and evolution relations. 

4.2.1 Elastoplastic kinematic decomposition 

The deformation at a material poinlQ is a combination of the elastic and plastic parts. In addition, 
a rotation occurs as part of the complete motion. These are shown schematically in Figure The 

^In the context of FEpX, a material point is a volume of material that is small in comparison to the finite element 
in which it resides (and thus, small in comparison to an individual crystal), yet large enough to fully reflect the crystal 
structure and the deformation processes. For slip, this means that the dimensions of the volume are much larger than a 
Burger’s vector. 
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Figure 2: Kinematic decomposition for motion by a combination of plastic slip, rotation and elastic 
straining. 


decomposition that describes this motion consists of breaking the deformation gradient into a sequence 
of three parts: a plastic part, a rotation and an elastic part, given as: 

f = f^rfP = v^r*fP (10) 

Each part of the decomposition brings the material point to a new configuration, starting with reference 
coordinates, X, and finishing at the current coordinates, x. The elastic part is a pure stretch which, 
by assuming small elastic strains, can be approximated with: 

v^ = l + e^ ( 11 ) 

where 

||e^||«l (12) 

The plastic part involves both stretch and rotation as a consequence of being a linear combination 
of slip modes, each of which is simple shear. The distinct rotation in f* (or equivalently, r*) is the 
rotation beyond that included in fP that is needed for consistency with the overall mapping given by 
Equation [T} 

The primitive solution field variable of FEpX is the velocity, owing principally to the code’s legacy 
of being a tool for modeling viscoplastic flow. To cast the kinematic decomposition in rate form, the 
velocity gradient first is written in terms of the deformation gradient and its time-rate-of-change: 

/ = ff-^ (13) 

where the velocity gradient is subsequently decomposed into the deformation rate and spin, as per 
Equation The deviatoric deformation rate is obtained by subtracting the volumetric part from the 
total: 

d' = d — d (14) 

Substituting Equation |10| into Equation and separating the deformation rate into its volumetric and 

deviatioric parts with Equation |14| gives: 


and 


tr(c/) = tr(e®) 


(15) 


d' = e^' + dP + e^'wP - wPe^' 


(16) 
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in which the small elastic strain approximation from Equation has been invoked. The spin becomes: 

w = + e^' dP — dP e^' (17) 

where, again, small elastic strains are assumed. In Equations and [I^ the hat over w/P and dP' 
indicates mapping to configuration B using r*, as indicated in Figure 

wP = r*wPr*^ (18) 

c/V = r*dP'r*^ (19) 

The deformations associated with elastic and plastic parts of the kinematic decomposition are inti¬ 
mately connected to the crystallographic lattice. The orientation of the lattice relative to a set of 
global base vectors at a designated material point is given by the associated quaternion, q, which is 
parameterized by the components: {qo,^- It is convenient to use other representations for orientation 
as well, depending on the task at hand. Two frequently used representations are the Rodrigues vector: 


and the rotation tensor: 


r = ntan | = q/qo 


R = 


1 + r ■ r 


(^1 {1 — r ■ r) + 2(r (g) r -|- / x r)) 


( 20 ) 


( 21 ) 


Properties that are dependent on lattice orientation are indicated as a function of q. Note also that 
changes in the lattice orientation that accompany the motion are embedded in r*, so the evolution rate 
of q is dehned in terms of the rate of change of r*. With the kinematic decomposition specified, the 
Kirchhoff stress is written based on the material point volume in the B configuration as: 


T = jdcr where fd = det(v'^) 


( 22 ) 


4.2.2 Fixed state constitutive relations 

The kinematic decomposition must be accompanied by equations relating the elastic and plastic defor¬ 
mations to the stress. For the elastic deformations, this relation is simply Hooke’s law, written using 
the B configuration as a reference volume: 


T=C{q)e^ (23) 

Here, the anisotropic behavior stemming from the crystal symmetry is indicated by the orientation 
dependence of the elastic stiffness. The structure of the elastic stiffness (occurrence of zero or repeated 
components in C{q)) reflects the application of symmetry conditions to the fully anisotropic version 
of Hooke’s law. These are more easily presented using a vector representation of the stress and strain 
tensors, as is done in Section 

For the plastic flow, the relation is a combination of several equations that describe crystallographic 
slip on a limited number of slip systems (commonly called restricted slip). First are equations for the 
kinematic decomposition written in terms of slip using the Schmid tensor’s symmetric and skew parts: 

iP = dP + wP (24) 

where 

dP = 7 "p" and h/?* = ^ 7 “^“ (25) 

OL OL 
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and 


p" = p“(q) = sym (s“ (g) m") (26) 

= q"(q) =skw(s"®m“) (27) 

The slip systems commonly assumed for face-centered cubic and hexagonal close-packed crystal types 
are shown in Figures and respectively. 



Figure 3; Primary slip systems for face-centered cubic (FCC) and body-centered cubic (BCC) crystal 
types. 



Pyramidal 
m = {10il] 
s = <1123> 



Basal 
m = {0001} 
s = <1120> 



Prismatic 
m = {10i0} 
s = <1120> 


Figure 4: Primary slip systems for hexagonal close-packed (HCP) crystal types. 


Next is an equation that defines the kinetics of slip, which introduces the rate dependence of plastic 
flow using a power law expression between the resolved shear stress and the slip system shearing rate: 


7" = /(T“,g“) = 70 


sgn(r° 


(28) 
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Crystal Type 

Name 

Number 

m 

s 

FCC 

octahedral 

12 

{111} 

< no > 

BCC 

- 

12 

{110} 

< 111 > 

HCP 

basal 

6 

{0001} 

< 1120 > 

- 

prismatic 

6 

{1010} 

< 1120 > 

- 

pyramidal 

6 

{1011} 

< 1123 > 


Table 1: Slip system vectors for (unstressed) FCC, BCC and HCP crystal types given a coordinate 
system attached to the lattice orientation. 


The resolved shear stress is scaled by the slip system strength, g°, which in general may be different 
for each slip system. However, in FEpX the slip system strengths are the same for each family of 
slip systems within a grain. Thus, the slip system strengths are all the same within each of the finite 
elements that discretize a grain for either FCC and BCC crystal types. For HCP crystals, the basal, 
prismatic and pyramidal strengths within each of the finite elements discretizing a grain can have 
different values, but all the systems of a given type have the same value. The values of the strength 
evolve with deformation according to the evolution equations, as discussed in Section 4.2.3 The ‘sgn’ 
term forces the shearing to be in the same direction as the shear stress. Finally, the resolved shear 
stress is the projection of the crystal stress tensor onto the slip plane and into the slip direction, which 
is readily computed with the symmetric part of the Schmid tensor: 


r“ = tr(p“T') 

The equations for slip are combined in a single, nonlinear relation as: 

dP =A4(q, 7 ") t ' 


(29) 


(30) 


Combining Equations 15, 16 


^ and into a single equation that relates the Cauchy stress to the 

This step is 


total deformation rate requires an additional step to discretize the elastic strain rate, 
introduced in Section [H 


4.2.3 State evolution equations 

There are two state variables at every material point that are updated as a deformation progresses, 
the lattice orientation and the slip system strength (also called hardnesses)]^ The rate of lattice re¬ 
orientation follows directly from the Equations |24| and |25[ assuming that the slip system shearing rates 
are known. Written in terms of the Rodrigues vector: 


r = -u> + (u> 
2 ^ 


r)r + (jj X r 


where 


cj = vect I ^ 


(31) 


(32) 


Evolution of the slip system strengths is governed by an additional, empirical relationship which 
follows as modified Voce form: 


5 " = ho 


gsii) - go 


7 


(33) 


^One could argue that the elastic strain also is a state variable, as it quantifies the shape of the unit cell. However, it 
is updated as an integral part of solving for the velocity field, rather than separate from solving for the velocity field as 
are summarized in this section for the lattice orientations and slip system hardnesses. See Section 
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Here, a saturation strength appears that is assumed to depend on a net local plastic strain rate 
computed from the sum of the magnitudes of the slip system shearing rate: 

( • \ m! 

and 7 = X] 

This equation is used to update the strength of the slip systems, family-by-family, in each element of 
the finite element mesh used to discretize a polycrystal. 


4.3 Boundary conditions 


Consistent with solid mechanics theoretical framework, the boundary condition applied to a surface of 
a virtual polycrystal may be either an imposed velocity or an imposed traction. For tractions, this is 
stated simply as; 

t{x) = t (35) 


while for the velocity condition, it is; 


v{x) = V 


(36) 


where the overbar indicates a known quantity. 


17 



5 Finite Element Implementation 


5.1 Matrix notation for tensorial quantities 


To facilitate the presentation and implementation of the finite element formulation, tensor quantities 
are written as matrices. Vectors map directly to one-dimensional column or row matrices. For second 
order tensors, column vectors are defined for with a designated ordering to the components. For the 
Kirchhoff stress and elastic strain, which are symmetric tensors, we use: 


r {r} = |rii T22 m V2T23 V2 tsi 
e" ^ {e®} = 6^2 6^3 \/2e|3 \/ 2 e^i \/2ef2} 


(37) 


(38) 


where the \/2 factor appears for the shear components in both tensors, which preserves the inner 
product relation 


T • = {rV {e^} 


(39) 


For the deviatoric parts of the second order tensors, a five-component form is adopted. For the Kirchhoff 
stress and the deformation rate; 


— I ”^ 22 ) v ^'^12 


(40) 


d' 


- I ^22) 


1 


— j ~ ^22) 


where the inner product again is preserved: 


^d '33 V2dU \/ 2 d' 3 i V2d 


-e' 

2^33 


^23 


31 


U 2 


V^^23 V^®31 V^®12 


■ d'= {r'f 


(41) 

(42) 

(43) 


Fourth-order tensors, namely the crystal elastic stiffness and compliance tensors, are commonly written 
as 6x6 matrices and populated according to the crystal symmetries. Hooke’s law written using the 
matrix form in a crystal coordinate bases for cubic and hexagonal crystal types are: 


and 


Til 


■ Cii C 12 C 12 


eii 

T 22 


C 12 Cn C 12 


622 

T 33 

V - 

C 12 C 12 Cn 

j 

633 

\ 

T 23 

' - 

C 44 

\ 

2623 

T 13 


C 44 


2613 

, n2 . 


C 44 


, 2ei2 , 


Cubic Symmetry 


(44) 


Til 


■ Cn C12 ^13 


611 

T22 


C12 Cn Ci 3 


622 

T33 

V - 

Ci 3 Ci 3 C33 

j 

633 

T23 

' - 

C44 

\ 

2623 

T13 


C44 


2613 

, T12 , 


(Cn - Ci2)/2 


, 2612 . 


Hexagonal Symmetry (45) 
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where only the nonzero values are displayed. A cautionary remark is added here regarding a factor 
of 2 that may appear with (744 in other expressions of Hooke’s law. FEpX expects values for the 


elastic moduli consistent with Equation 44 or even though stiffness or compliance matrices internal 
to FEpX are constructed somewhat differently. A restriction is placed on the moduli for hexagonal 
crystals vis-a-vis coupling of the shear and volumetric responses, as described in the next paragraph. 

In the kinematic development, the motion is split into volumetric and deviatoric parts according 
to Equations 15 and 16 This split is convenient for the numerical implementation, but limits the 


generality of the hexagonal behaviors. In particular, the volumetric and deviatoric responses separate 
only when Cn + C 12 = Cia -|- 6*33. Thus, only four of the five nonzero moduli are independent. To 
guarantee that this constraint is imposed, (Cn, C 12 , C13, C44) are read in the input for FEpX and C33 
is computed to satisfy the constraint. 

When split into volumetric and deviatoric parts. Equation gives: 


and 


tr{r} = -tr{e'’} 


{r'}= c' {e-} 


(46) 


(47) 


where the stress and strain vectors are consistent with Equations 37 - Table lists the values of 


the K and the nonzero (diagonal) entries of d in terms of the moduli presented in Equations 


44 


and 


Parameter 

Cubic 

Hexagonal 

K 

3(C'ii + 2 C 12 ) 

3(Cii + C12 + C13) 

Cn 

Cii - C12 

Cii - C 12 

^22 

Cii - C 12 

3(^33 - C13) 

C33 

C44 

(744 

C44 

C44 

C44 

*^55 

(744 

Cii - (7 i2 


Table 2: Values of the moduli used in the separated form of Hooke’s law 


5.2 Time-discretized elastoplastic relations 


Equations pAj [Tbl andare now merged into a single equation that relates the Cauchy stress to the 
total deformation rate. First, the spatial time-rate change of the elastic strain is approximated with a 
finite difference expression: 

{e“} = ^ ({<='}-{«) (48) 

where {e®} is the elastic strain at the end of the time step and {ep} is the elastic strain at the beginning 
of the time step. The difference approximation is employed in an implicit algorithm, wherein the 
equations are solved at the time corresponding to the end of the time step. This time corresponds to 
the current configuration. Writing the time rate change of the strain in terms of strains at two times 
facilitates substitution of Hooke’s law - namely at the end of the time step. The elastic strain at the 
beginning of the time step is known from the solution for the preceding time step. For the volumetric 


part of the motion, combining Equations 15 and 23 with the difference expression gives: 




— vr = 


/3 


tr{d} + ^tr{e^} 


(49) 
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Turning to the deviatoric (shearing) part of the motion, inserting Equation 48 into Equation 16 gives: 

{d'l = T {eq + {d^} + 




(50) 


where 




is the matrix form of w^: 



■ 0 

0 

-2^2 

-<3 

^23 



0 

0 

0 

y/SujP^ 

\/3W2^ 

wP 

= 

2^2 

0 

0 

-^23 




^13 



0 

-WP^ 


- -^23 

— \/3h)23 

''P 

<3 

^p 

<2 

0 


(51) 


The equations for plastic slip (Equation |30[) for the plastic deformation rate: 


{d^} 


m 


m 


E 




(52) 

(53) 


together with Equation]^ are substituted to render an equation that gives the deviatoric Cauchy stress 
in terms of the total deviatoric deformation rate and a matrix, |h|, that accounts for the spin and the 
elastic strain at the beginning of the time step: 


P'} = [s]({d'}-{h}) 


where: 


1 -1 


A 

At . 


1 -1 


+ /3 


m 


{^} 






(54) 

(55) 

(56) 


Equations 54 55 and 56 will be used in the weak form of equilibrium to write the stress in terms of 
the deformation rate. 

5.3 Interpolation functions 

FEpXemploys a standard isoparametric mapping framework for discretizing the problem domain and 
for representing the solution variables. The mapping of the coordinates of points provided by the 
elemental interpolation functions, N(^,t/, <^) , and the coordinates of the nodal points, {X}: 


{x}= N(e,r?,C) {X} 


(57) 


where (^, rj, C) are local coordinates within an element. The same mapping functions are used for the 
solution (trial) functions which, together with the nodal point values of the velocity, |v|, specify the 
velocity held over the elemental domains: 


{u}= [n(C,7?,C)]{v} 


(58) 
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The deformation rate is computed from the spatial derivatives (derivatives with respect to x) of the 
mapping functions and the nodal velocities as: 




(59) 


B 


is computed using the derivatives of 


with respect to local coordinates, to- 


N(^,f?,C) 

gether with the Jacobian of the mapping specified by Equation following standard finite element 
procedures for isoparametric elements. 

FEpX relies principally on a 10-node, tetrahedral, serendipity element, as shown in FigureThis 
element provides pure quadratic interpolation of the velocity field. FEpX employs a Galerkin 





Figure 5: 10-node tetrahedral element with quadratic interpolation of the velocity, shown in the parent 
configuration and bounded by a unit cube. 


methodology for constructing a weighted residual. The weight functions therefore use the same inter¬ 
polation functions as used for the coordinate map and the trial functions: 


1^^} = [N(e,r?,C)]{iZ'} 


(60) 


5.4 Finite element residual for the velocity field 

Equilibrium is enforced by requiring a global weighted residual to vanish: 

Ru = I V' • (divcr^ -|- t) dB = 0 (61) 

Jb 

The residual is manipulated in the customary manner (integration by parts and application of the 
divergence theorem) to obtain the weak form: 

Ru = — tr gradi/’) dB + tt div0d;B + t ■ ■j/’dT + l ■ ijjdB (62) 

Jb ^ ^ Jb Jsb Jb 
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Introduction of the trial and weight functions gives a residual vector for the discretized weak form for 
each element: 


{Ru^} 


where 


i^ele 


^ele 
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,ele 


{v} - {ff} - - {f*} 
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iT /3 
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f [n(c,7?,c)]^{4 

{<l}^{eS}d^ 


dB 


{f*} 


B 


1 T r 1 T r 1 

X 


{^} 


di3 


(63) 

(64) 

(65) 

( 66 ) 

(67) 

( 68 ) 


The integrals appearing in Equations 64 68 are evaluated by numerical quadrature. 

Assembling the elemental matrices and requiring that the residual vanish for all independent vari¬ 
ations in the weights gives: 




+ 


K,. 


{v} = {f.} + {f.} + {f.} 


The ess enti al boundary conditions are applied prior to 


Section 


5.6 


The matrix equation given by Equation 


69 


solving for the nodal ve 


is nonlinear 


Kd 


and 


The solution methodology used in FEpX is outlined in Section 5.5 


(69) 

ocities, as described in 
depend on Iv|). 


5.5 Nonlinear solution algorithm for obtaining the velocity field 


To solve Equation 69 for the velocity field, an iterative methodology is invoked. This methodology 
is a hybrid procedure that utilizes a combination of successive approximations (Picard) and Newton- 
Raphson updates. To accomplish this, the assembled residual force vector, |i?u|, is defined as: 




K, 


+ 


K, 


{v}-{f„}-{f,}-{f„} 


(70) 


The goal of the iterative process is to drive the residual to zero through a series of corrections, |au|, 
to an estimate of the velocity field. Denoting the estimate of the velocity on the iteration as |v| 
and the estimate on the next iteration as 1 V > , the iteration procedure is written simply as: 


V 




i+1 


where 


r I 

{au} 


is determined from the solution of 

type 


K. 


+ 


K. 




Here, 


K. 


1 type 


refers to either a tangent modulus. 


K. 


-| tan 


or a secant modulus. 


K. 


(71) 


(72) 


, as specified 


by the hybrid procedure. Convergence is based on changes in the norm of |v| becoming small. 
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5.6 Time marching and bonndary conditions 


The intent of the FEpX framework and the derivative finite element code is to simulate the deformation 
of virtual polycrystals over time. To this end, time histories are approximated by solving for the velocity 
field at a series of discrete times. Simply stated, FEpX computes the velocity field at the end of time 
interval using Equation 70 knowing the velocity field and state at beginning of the time interval. The 
geometry and state variables (lattice orientations and slip system strengths) are updated concurrently 
with the velocity field at the end of the time step. The time marching method is documented in mm- 

Over the course of a deformation history, the applied boundary conditions often change. This may 
imply that either natural or essential boundary conditions in Equations |35| and |36| are functions of time. 
Presently, FEpX allows the user to change the values of imposed velocities or forces at nodes over the 
course of a deformation, but does not allow the user to change the type of boundary condition. That 
is, a nodal point with imposed velocity (essential boundary condition) will have an essential boundary 
condition throughout a simulation, but the value of the velocity that is imposed may change with time. 
The same applies for nodes with natural boundary conditions - the force may change with time, but 
the condition at that node will remain a natural boundary condition throughout the simulation. 

Given this limitation, the boundary condition options within FEpX accommodate a number of 
possibilities. For example, it is anticipated that a common use of FEpX will be to simulate the response 
of virtual polycrystals being subjected to boundary conditions that replicate mechanical tests performed 
using a load frame. The control of mechanical tests can be designed to provide programmed force 
histories, programmed displacement histories, or combinations of these. FEpX offers the capabilities 
to impose boundary conditions to mimic mechanical test histories. Some of the possibilities include; 


• specified crosshead/actuator velocity, 

• specified load history, 

• unloading episodes, and 

• uniaxial and biaxial loading modes. 
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6 Input and Output Data 


The modeling framework described in this article is intended for simulating the deformations of poly¬ 
crystalline aggregates. The material state (including phase topology, grain morphology, crystallographic 
texture, slip system strengths) both influences the deformation and is affected by the deformation, and 
thus is an integral part of the definition of a material system being modeled. The finite element for¬ 
mulation provides the solution methodology for solving the system of equations under imposed loading 
for the polycrystalline aggregate of interest. Thus, the input and output data for FEpX are similar 
in content to other finite element codes, but are tailored for and organized around the polycrystals. 
One must define the finite element mesh, assign material attributes consistent with the state to the 
elements, and specify boundary and initial conditions. In addition, information needed for internal 
computational procedures, such as choices regarding type of nonlinear solver to employ, convergence 
tolerances on nonlinear iterations, as well as the frequency and extent of output data must be provided. 
A general description of the input data is given here. The detailed description of the input data needed 
by a user to execute FEpX is given in a separate document. 

Because FEpX uses microstructural information generated by experiment or simulation to assign 
material attributes, pre-processing routines are needed to prepare input files. Other software package 
serve this function. OdfPf is a suite of Matlab routines written to perform a variety of tasks related 
crystallographic texture and is recommended as a complement to FEpX. Neper is capable of meshing 
virtual polycrystals created by Voronoi tesselation and constructs the nodal point and element data 
needed by FEpX. Instructions for using these packages are included with the respective packages. 
Other packages for visualization, interpretation and archiving data also are needed for post-processing 
FEpX results. Again, OdfPf provides capabilities to translate FEpX output files into formats 
appropriate for using other software packages. 

6.1 FEpX input data 

The body being loaded and deformed in a FEpX simulation is a virtual polycrystal - a set of grains 
(each grain being a single crystal) that forms a fully-dense solid. The finite element mesh needed to 
build a virtual polycrystal must be created in advance, which can be done with Neper or other mesh 
generation codes. The input data is organized as follows: 

• Setting up a job — a script to execute FEpX together with a file which lists the input files 
described below. 

• Defining a virtual polycrystal — a set of files: a file that specifies the single crystal elastic 
and plastic material properties for each phase; a file that defines the finite element mesh (nodal 
coordinates, element connectivity, and surface elements); a file that designates the phase and 
grain numbers for each element; a file that provides the starting lattice orientation for each grain; 
and, one or more files that define the vertices of the single crystal yield surface for each phase. 

• Controlling the deformation — a set of two files: a file designates the type of boundary 
condition (essential or natural) for each degree of freedom of all the nodes and the corresponding 
velocity or force; and, a file specifies target loads or displacements, depending on the loading 
type. 

• Postprocessing for diffraction data — a file that provides information for averaging strains 
and stresses over crystallographic fibers. 

• Prescribing optional input — a file that specifies the choices for loading protocols, post¬ 
processing options, various solution method options, and associated convergence limits. 
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6.2 FEpX output data 

In the present implementation, the FEpX code writes files for solution variables at times designated 
in the input data. The solution variables written to files are: 

• Geometry: the nodal point coordinates and the nodal point velocities. 

• Deformation: the effective deformation rate, the effective plastic deformation rate, the slip 
system shearing rates, the effective strain, and the effective plastic strain. 

• State Variables: the slip system strengths and the lattice orientations. 

• Stress: the stress tensor. Note that the components of the Cauchy stress tensor are written in 
the global frame in the following order: uii, (T12, uis, 1T22, (^ 23 , 0 - 33 . 

• Elastic strain: the elastic strain tensor. As with the stress, the components of the elastic strain 
tensor are written in the global frame in the following order: en, ei2, ei3, 622, £23, €33- 

• Fiber averages: mean and standard deviation values of the stress, elastic strain, and slip system 
activities taken over designated crystallographic fibers. 

• Monitors: the resultant force on each external surface. 

When executed on a parallel architecture, FEpX writes output files on every node which are written 
to the master node at the end of the job. This leaves the results spread over many output files and 
in a format that is inconvenient for postprocessing. An OdfPf script is available to concatenate all of 
the output files into a single data structure that is readily used within postprocessing tools. 

6.3 Exporting Input and Output Data 

Exporting input and output data is done for two reasons. One is so that the results may be archived 
in a data management system appropriate for the project of interest. The other is to interface with 
visualization codes or other specialty post-processing software, such as forward projectors or virtual 
instruments. The OdfPf package includes scripts for these purposes: 

• The HDF-5 file format is a commonly used standard for writing archivable files. A matlab script 
is available that prepares an HDF-5 file with the FEpX standard input and output data. This 
script can be modified to include other information (such as postprocessing results) if desired. 

• To facilitate visualization of the results, matlab scripts are available to export files that can be 
read in by Paraview or Visit. The applications of FEpX shown in Section were plotted with 
Paraview. 

6.4 OdfPf capabilities 

OdfPf figures prominently in the use of FEpX to simulate the behavior of virtual polycrystals, both 
in the instantiation of virtual polycrystals and in the comparison of simulation results to diffraction 
data. OdfPf is a function set is a collection of MATLAB functions which operate on ODF’s (orien¬ 
tation distribution functions) and PF’s (pole figures). It handles plotting of the ODF using Rodrigues 
parameters, plotting of pole figures and inverse pole figures, evaluation of pole figures inverse pole 
figures from ODF’s, and it provides many tools for computing ODF’s from pole figures. Archival pub¬ 
lications are available that cover various aspects of its use or the use of a Rodrigues parameterization 
of orientation space in quantitative texture analyses. Relevant articles include: 
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7 Example Problems 


To illustrate the use of FEpX in simulating the mechanical response of virtual polycrystals, three 
example problems are presented. We use different methods for instantiating virtual polycrystals for 
each of the three examples. The first is an example of building a virtual polycrystal comprised of 
regular, dodecahedral-shaped grains. In the second example, the virtual polycrystal is a discretized 
Voronoi tessellation. The third example demonstrates the construction of a virtual polycrystal by 
mapping serial sectioning data directly onto a uniform mesh. All examples were executed on a parallel 
computer using 64 processors (8 nodes each with 8 cores). The same version of FEpX was used for 
all examples with data files conforming to the input manual. Postprocess of the results was performed 
with OdfPf scripts to write *.vtk files. The *.vt files were imported into Paraview for plotting. 

7.1 Virtual Polycrystal Generated by Regular Tessellation 

This example was provided by Andrew Poshadel and was developed as part of his research on yielding 
of dual phase alloys [l 6 ] . The example demonstrates the application of FEpX to a virtual polycrystal 
with two phases that was built using dodecahedral grains. The two phases are the austenitic and 
ferritic phases of a dual steel (LDX-2101), referred to symbolically as the 7 and a phases, respectively. 
The stock material has a microstructure consistent with having been rolled or extruded. 


7.1.1 Defining the virtnal polycrystal 


Following the summary of the input data given in Section 6.1 
FEpX can be organized into five categories: 


the required input data to execute 


1. Phase Attributes: The a and /3 phases both have cubic crystal structure - FCC for the 7 phase 
and BCC for the a phase. The single crystal, cubic, elastic moduli for the two phases are listed 
in Table Based on experimental data that indicates the plastic behaviors of the two phases 
are comparable, the same plasticity parameters were assigned to both phases for the purpose of 
this example. These parameters are listed in Table The slip systems are different for the two, 
however, with the FCC 7 phase using the {111} < 110 > systems and the BCC a phase using 
the {110} < 111 > systems (See Table[^. This information is provided to FEpX in the *.matl 
input file. 


Phase 

Type 

Cii (GPa) 

C12 (GPa) 

Cm (GPa) 

7 

FCC 

204.6 

137.7 

126.2 

a 

BCC 

236.9 

140.3 

116.0 


Table 3: LDX 2101 elastic moduli used in the single crystal constitutive equations for the two-phase 


virtual polycrystal. Values listed conform to the convention given in Equation 44 


Phase 

0 

1 

m 

ho (MPa) 

go (MPa) 

n' 

91 (MPa) 

Is (s ^) 

m' 

7 

1.0 

0.02 

391.9 

237.0 

1 

335.0 

5.0 X 10^^ 

0 

a 

1.0 

0.02 

391.9 

237.0 

1 

335.0 

5.0 X 10^° 

0 


Table 4: LDX 2101 slip parameters used in the single crystal constitutive equations for the two-phase 
virtual polycrystal. Values listed conform to the convention given in Equations [Mj and 34 
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2. Mesh definition: A finite element mesh underlying the virtual polycrystal was instantiated 
using a custom MATLAB script. It consists of a regular, rectangular layout of elements. These 
elements can be grouped to form regular dodecahedral grains. The resulting mesh, shown in 
Figure has 117,504 10-node tetrahedral elements and 173,829 nodal points. The corresponding 
arrays for the nodal point coordinates and element connectivities are provided to FEpX in the 
*.mesh input file. 



Figure 6 : Finite element mesh for the dual phase steel virtual sample. 


3. Phase and grain definition: The finite elements must be assigned to phases and grains. 
The first step was to define the spatial distributions of the two phases in one plane using a 
regular layout of dodecahedral grains. The second step was to extrude (repeat) the planar 
layout in the direction perpendicular to the plane to create the cube-shaped polycrystal with 
a microstructure similar to the stock material. In this example, elements are one of the two 
phases (7 or a). Within subdomains of a single phase, there may be one or more crystals. For 
this dual phase steel, the subdomains of both phases have multiple crystals, although the 7 phase 
subdomains generally had fewer crystals than the subdomains of the /? phase subdomains. The 
grain definitions follow from the assignment of lattice orientations to the elements. Contiguous 
elements with the same orientation constitute a grain. Lattice orientations were chosen randomly 
from measured orientation distributions and assigned to element to create the desired grain 
arrangement within the phases. Each grain also was assigned the same initial slip system strength, 
which in turn was assigned to every element of the grain. Figure shows the phase and grain 
assignments associated with the mesh shown in Figure This information is provided in the 
*.grains and *.kocks files. 

4. Vertex files: standard definitions of the single crystal yield surface vertices were used for both 
the FCC and BCC phases. 
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(a) Phases (b) Grains 


Figure 7: Phase and grain distributions for the dual phase steel virtual sample. For the phase distri¬ 
bution, blue indicates the BCC ferritic (a) phase; red indicates the FCC austenitic (7) phase. For the 
grain distribution, grains are indicated by domains of uniform color. 

7.1.2 Controlling the loading: 

1. Boundary conditions: The boundary conditions are intended to simulate the loading applied 
in a tension test. The virtual polycrystal is constrained on the bottom and stretched in the z 
direction by an imposed axial velocity on the top. Two adjoining lateral surfaces are traction-free 
while a symmetry condition is applied on the other two. Rigid body translations and rotations 
have been suppressed by the application of the symmetry conditions. This information is given 
in the *.bcs file. 

2. Target loads: Simple load control is applied to extend the sample. Several z-direction target 
loads along a monotonically increasing path (no unloading episodes) are specified to provide 
points for writing output data. This information is given in the *.loads file. 

7.1.3 Specifying options: 

1. Load controls: The “control by load” mode is used to control the loading history. 

2. Convergence criteria: The default parameters have been used for both the velocity field and 
crystal stress iterative procedures. The Newton-Raphson procedure is invoked for the velocity 
solutions. 

7.1.4 Selected Simulation Results 

Simulation results are available for postprocessing at the points in the loading designated by the target 
loads. These results were aggregated and written to a .vtk file for plotting with Paraview. Figure 
shows the distribution of axial component of the stress and the effective plastic strain at the last target 
load of 590 N. At this load, the nominal axial strains was approximately 10%. The stress shows spatial 
variations due to the anisotropy of the crystal properties and the interactions among the grains. There 
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is an increase in the variability of the axial component of the stress as the stress level is increased due 
in part to the change in the principal directions of the crystal stresses as the stresses move toward a 
vertex of the single crystal yield surface during the elastic-plastic transition. The effective stress (not 
show here) exhibits less variability. Figure shows the evolution of strength over the course of the 
loading. The deformation of the mesh has been exaggerated by a factor of 4 to facilitate visualizing 
the heterogeneity of the deformation. A correlation between the plastic deformation and the strain 
hardening is evident. A great deal more information is available in the simulation output. For 




(a) Axial stress component. (b) Effective plastic strain. 

Figure 8: Axial stress and plastic strain distributions at nominal load of 590N shown on the deformed 
mesh. 



Figure 9: Slip system strength at 590N shown on a exaggerated (x4) deformation field. 
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example, lattice strains in crystals on designated crystallographic fibers was collected and averaged for 
comparisons to experiments in which lattice strains were measured by neutron diffraction during in 
situ loading. 


7.2 Voronoi Tessellated Virtual Polycrystal Generated with Neper 


This example was provided by Amanda Mitch and was developed as part of her research on reduced- 
order representation of crystal stress distributions for use in a methodology for quantifying residual 
stress distributions in engineering components m- The example demonstrates the application of 
FEpX to a virtual polycrystal that was built using a Voronoi tessellation to define the grains. The 
material is single phase and the single crystal properties are consistent with a FCC crystal type, but 
do not represent any particular metal or alloy. The input data, following the organization given in 
Section 16.11 is summarized below. 


7.2.1 Defining the virtnal polycrystal 

1. Phase Attribntes: The material is single phase with a cubic (FCC) crystal structure. The 
single crystal, cubic, elastic moduli are listed in Table Plasticity parameters are generic, 
being similar to copper alloy. These parameters are listed in Table The slip systems are the 
customary primary systems for FCC crystals: the {111} < 110 > systems. This information is 
provided to FEpX in the *.matl input file. 


Phase 

Type 

Cii (GPa) 

Ci 2 (GPa) 

Cm (GPa) 

a 

FCC 

245. 

155. 

62.5 


Table 5: Elastic moduli used in the single crystal constitutive equations for the Voronoi-based virtual 
polycrystal. 


Phase 

O 

1 

m 

ho (MPa) 

go (MPa) 

n' 

91 (MPa) 

is (s ^) 

m' 

a 

1.0 

0.05 

200. 

210. 

1 

330. 

5.0 X 10^*^ 

5.0 X 10“^ 


Table 6: Slip parameters used in the single crystal constitutive equations for the Voronoi-based virtual 
polycrystal. 

2. Mesh definition: A virtual polycrystal was instantiated using the Neper code. Neper builds 
a Voronoi construction of the domain to define grains and then discretizes the grains into hnite 
elements. The resulting mesh, shown in Figure [Tol has 96,758 10-node tetrahedral elements and 
134,362 nodal points. The corresponding arrays for the nodal point coordinates and element 
connectivities are provided to FEpX in the *.mesh input file. 

3. Phase and grain definition: All grains are the same phase. The grain designations for the 
finite elements follow from the Voronoi tessellation and are shown in Figure [T^ This information 
is provided in the *.grain and *.kocks files. 

4. Vertex files: a standard definition for the vertices of a FCC single crystal yield surface was 
used. 
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(a) Mesh (b) Grains 

Figure 10: Mesh and grains for the iVeper-built polycrystal. 


7.2.2 Controlling the loading: 

1. Boundary conditions: The boundary conditions are intended to simulate the triaxial loading 
of the sample such that the stress components remain in constant proportions given by.: 


[o-] = 0-1 


1.0 0 0 

0 -0.625 0 

0 0 -0.375 


(73) 


Normal velocity components of different magnitudes are applied to three adjacent surfaces while 
the opposing surfaces are hxed in place. The surface velocities are adjusted to achieve tractions 
consistent with the target ratios for triaxial stress state. 

2. Target loads: Two target loads along a monotonically increasing path (no unloading episodes) 
are specified to provide points for writing output data. At the first target load, the response is 
essentially elastic; the second target load is sufficient to induce plastic strains on the order of 3%. 
The target load information is given in the *.loads file. For triaxial loading, three normal forces 
are specified for each target load consistent with the desired stress state specified in Equation 

7.2.3 Post-Processing: 

1. Lightup: Fiber-based quantities are computed for 6 fibers: (100) and (111) crystal planes in the 
[100], [010]and[001] sample directions. 


7.2.4 Specifying options: 

1. Load controls: The “Triax-CSR” mode is used to control the loading history. Using this option, 
the magnitude of the velocities are adjusted to impose a specihed loading rate (increase in the 
applied forces). The relative values of the imposed surface velocities are adjusted to maintain the 
specified state of triaxial stress. 
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2. Convergence criteria: The default parameters have been used for both the velocity held and 
crystal stress iterative procedures. The Newton-Raphson procedure is invoked for the velocity 
solutions. 

3. Lightnp: The option to compute hber-based quantities is specihed. 

7.2.5 Selected Simulation Results 

Stress distributions for the polycrystal are shown in Figure [TT| at the second target load. In these images, 
the normal components of the stress are plotted over the deformed mesh. The total deformation is 
not large, so the change in shape from the initial conhguration shown in Figure is difficult to 
discern. The differences in overall shade between the three images reflects the triaxial stress condition 
intentionally imposed on the polycrystal. There are variations over the polycrystal for all of the 
components stemming from the elastic and plastic anisotropy. 



(a) XX component (b) yy component 



(c) zz component 


Figure 11: Normal stress component distributions at an = 


In Figure 12 the normal components of the elastic strain are depicted. 


225MPa 

A noticeable contrast to 
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the stress distributions is evident. For the elastic strains, the net effect of the stress triaxiality and 
grain interactions is to produce distributions that span approximately the same range in strain for all 
the normal strain components. Unlike the stress, it is difficult to discern the triaxiality of the stress 
from differences in the lattice lattice (elastic) strain distributions. Figure [I^ shows the effective plastic 
strain at the second target load. Again, the distribution is not uniform: some elements display several 
percentage plastic strain while other have almost no plastic strain. Details are available in m- 



(a) XX component 


(b) yy component 



(c) zz component 


Figure 12: Normal lattice (elastic) strain component distributions at ai = 225MPa 

The distribution of stresses over a polycrystal depends strongly on the levels of anisotropy in the 
elastic and plastic behaviors of the constituent crystals. The stress state must vary spatially to satisfy 
compatibility and equilibrium if the properties vary. Often these level are different, which is evident 
by observing the stress distributions as the polycrystal is loaded through the elastic-plastic transition. 
At lower loads, the behavior is essentially purely elastic and the distribution is controlled by the elastic 
moduli. At high loads, the behavior is dominated by the single crystal yield surface. This is illustrated 
by the relative changes in average elastic strains along the selected fibers given in Table One can 
readily observe that between the target loads for ai = 200MPa and ai = 225MPa there is not the same 
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Figure 13: Effective plastic strain distribution at ai = 225MPa 


proportional increase in strains, as would be expected in the response remained linear and superposition 
could be applied. The adjustment in stress that accompanies yielding implies that the principal axes 
of the strain rotate as the stress moves toward a vertex of the single crystal yield surface. 


cTi (MPa) 

(100)||[100] 

(111)||[100] 

(100)||[010] 

(111)||[010] 

(100) II [001] 

(111)||[001] 

200 

0.00211 

0.00196 

-0.00115 

-0.00117 

-0.00075 

-0.00078 

225 

0.00221 

0.00244 

-0.00105 

-0.00160 

-0.00044 

-0.00113 


Table 7; Average lattice strains along selected fibers. 
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7.3 Voxel-Based Virtual Polycrystal Generated from 3-D Serial Section Maps 

This example was provided by Donald Boyce and was developed as part of an ONR-sponsored project 
on strength and ductility in titanium alloys. The example demonstrates the application of FEpX to 
a virtual polycrystal that was built by mapping voxel data to a regular mesh to define the grains. The 
titanium alloy being modeled is Ti-6A1-4V, which is two-phase at room temperature. However, since 
the volume fraction of the BCC phase is only about 7%, this analysis is performed assuming a single 
(HCP) phase. 

7.3.1 Defining the virtnal polycrystal 

1. Phase Attributes: The principal phase for this titanium alloy has hexagonal symmetry (HCP), 
and is designated as the a phase. The single crystal, hexagonal, elastic moduli for it are listed in 
Table Input to the code does not include C33. Rather, it is computed internally to assure that 
the constraint to decouple the deviatoric and volumetric responses is satisfied. The plasticity 
parameters were estimated from fitting stress-strain data for the alloy and are listed in Table 
This type of alloy exhibits very little strain hardening; the choice of parameters accomplishes 
this by setting the initial slip system strength to the saturation value. The slip systems include 
prismatic, basal and pyramidal systems as per Table This information is provided to FEpX 
in the *.matl input file. 


Phase 

Type 

Gii (GPa) 

C 12 (GPa) 

Gi 3 (GPa) 

Cm (GPa) 

a 

HCP 

161.4 

91.0 

69.5 

46.7 


Table 8; Titanium elastic moduli used in the single crystal constitutive equations for the voxel-based 
virtual polycrystal. 


Phase 

70 (s ^) 

m 

ho (MPa) 

go (MPa) 

n' 

91 (MPa) 

is (s ^) 

m' 

a 

1.0 

0.01 

1000. 

500. 

1 

500. 

5.0 X 10^^ 

0.01 


Table 9: Titanium slip parameters used in the single crystal constitutive equations for the voxel-based 
virtual polycrystal. 


Basal 

Prismatic 

Pyramidal 

1 

1 

3 


Table 10: Relative strength for the titanium slip system used in the single crystal constitutive equations 
for the voxel-based virtual polycrystal. 


2. Mesh definition: The finite element mesh underlying the virtual polycrystal was instantiated 
by mapping voxel-based (3D) orientation map onto a regular mesh using a custom MATLAB 
script (available in the OdfPf package). The orientation map was obtained from serial section 
data measured using electron back-scattered diffraction (EBSD). The finite element mesh spans 
a volume of 20pm x 20pm x 60pm that coincides with an interior portion of the experimental 
volume. The element size was chosen to give and resolution comparable to the spatial resolution 
of the data. The resulting mesh, shown in Figure [T4l has 144,000 10-node tetrahedral elements 
and 230,401 nodal points. The corresponding arrays for the nodal point coordinates and element 
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connectivities are provided to FEpX in the *.mesh input file, along with definition of the six 
sample surfaces in terms of the mesh elements. 



(a) Finite element mesh (b) Lattice orientations 


Figure 14: Finite element mesh for the mill annealed titanium alloy and grain orientations assigned to 
the elements. 

3. Phase and grain definition: For every element of the mesh, the voxel that is closest to the 
element is identified using the distance between the centroid of the element and the centroid of 
the voxel). The orientation data of the voxel is then used to assign the lattice orientation for 
the element. The simulation assumes all grains are of the same HCP phase. All the grains are 
assigned the same initial slip system strength, which in turn was assigned to every element of the 
grain. Figure [T4| also shows the grain assignment associated with the mesh. While some grains 
are evident, noisy or missing information in the orientation data makes crisp dehnition of grains 
difficult. The grain assignment and lattice orientation data are given in the *.grain and *.kocks 
files. 

4. Vertex files: Vertices of a HCP single crystal yield surface was used having a topology consistent 
with the prescribed ratios of the slip system strengths of (1 : 1 : 3) for the basal, prismatic and 
pyramidal slip systems, respectively. 

7.3.2 Controlling the loading: 

1. Boundary conditions: Boundary conditions are chosen to mimic a tensile test: there is a 
fixed velocity applied on the upper surface (positive y face) while the lower surface is hxed from 
translation in the y direction. The lateral surfaces have two traction free surfaces (positive x and 
z) and two symmetry surfaces (negative x and z). This information is in the *.bcs file. 

2. Target loads: Simple load control is applied to compress the sample. Three y-direction target 
loads along a monotonically increasing path (no unloading episodes) are specified to provide 
points for writing output data. The final target load was sufficient to compress the sample by 
approximately 1%, overall. This information is given in the *.loads file. 
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7.3.3 Specifying options (information in the *.options file): 

1. Load controls: The “control by load” mode is used to control the loading history. 

2. Convergence criteria: Default parameters have been used for both the velocity field and 
crystal stress iterative procedures. The Newton-Raphson procedure is invoked for the velocity 
solutions. 


7.3.4 Selected Simulation Results 


Figure 15 shows the axial stress and the axial lattice strain at the final target load. Here the grain 


structure is more evident. Because the load is sufficient to cause wide-spread yielding, higher stresses 
and strains typically occur in grains whose lattices are at stronger orientations. Note that the stress 
distribution is not merely a scaled version of the strain distribution. This is a result of the elastic 
anisotropy, which implies that the eigenvectors of the stress and strain tensors do not necessarily align. 

At first glance, the stress levels depicted in Figure [T5| might appear to exceed stress limits imposed 
the single crystal yield surface. However, the large grain size relative to the sample size has an effect 
such that the deformation is more highly constrained. The consequence is that the mean stress (shown 
in Figure [Th]) is larger at many points than would be expected for a case of simple tension. The plastic 
straining induced by this highly heterogeneous stress held is very inhomogeneous, as is evident from 
the distribution of the effective plastic strain shown in Figure 16 This plot shows how plastic flow 


interconnects through the polycrystal leaving some domains relatively undeformed. 




(a) Axial stress 


(b) Axial lattice (elastic) strain 


Figure 15: Axial stress and axial lattice strain distributions at the third target load. 
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(a) Mean stress (b) Effective plastic strain 

Figure 16: Mean stress and plastic strain distributions at the third target load. 
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